﻿using System;
using System.Numerics;

namespace WaveSim.Core.Equations.Galerkin
{
    public class Galerkin2DEquation : IWaveEquation, IEquation
    {
        private static Complex im = new Complex(0, 1);

        private WaveGuide waveGuide;
        private BaseDielectric dielectric;

        public Galerkin2DEquation(Wave wave, WaveGuide waveGuide, BaseDielectric dielectric, int nx)
        {
            Wave = wave;
            this.waveGuide = waveGuide;
            this.dielectric = dielectric != null ? dielectric : new Composite();
            this.Nx = nx;
        }

        public int Nx { get; set; }

        public WaveGuide WaveGuide
        {
            get
            {
                return waveGuide;
            }
        }

        public Wave Wave { get; set; }

        public Complex Get(int i, double x, Complex[] y)
        {
            if (i >= y.Length)
            {
                return default(Complex);
            }

            if (i % 2 == 0)
            {
                return y[i + 1];
            }
            else
            {
                int m = (i + 1) / 2;
                Complex Am = y[i - 1];
                double arg = m * Math.PI / waveGuide.A;
                Complex sum1 = arg * arg * Am;
                Complex sum2 = 0;
                for (int j = 0; j < y.Length; j += 2)
                {
                    Complex k0 = Wave.K(waveGuide, j / 2 + 1);
                    sum2 += y[j] * k0 * k0 * Integral(x, j / 2 + 1, m);
                }
                return (sum1 - sum2);
            }
        }

        private Complex Integral(double z, int n, int m)
        {
            double x0 = 0;
            double x1 = waveGuide.A;
            double hx = (x1 - x0) / Nx;

            //waveGuide.B;

            //The composite trapezoidal rule
            double arg = (x1 - x0) / Nx;
            Complex summa = (IntegralFunction(x0, z, n, m) + IntegralFunction(x1, z, n, m)) / 2;
            for (int i = 1; i < Nx - 1; i++)
            {
                double x = x0 + i * hx;
                summa += IntegralFunction(x, z, n, m);
            }
            return summa * hx * waveGuide.B;
        }

        private Complex IntegralFunction(double x, double z, int n, int m)
        {
            return dielectric.Eps(x, z) * dielectric.Mu(x, z) * Phi(x, n) * Phi(x, m);
        }

        private double Phi(double x, int n)
        {
            return Math.Sqrt(2 / (waveGuide.A * waveGuide.B)) * Math.Sin(n * Math.PI * x / waveGuide.A);
        }

        public Complex[] StartValues(double[] amps)
        {
            if (amps == null || amps.Length == 0)
            {
                return default(Complex[]);
            }

            int modeCount = amps.Length;
            var startValues = new Complex[2 * modeCount];
            for (int i = 0; i < modeCount; i++)
            {
                // i
                startValues[i] = amps[i];

                // i + 1
                int m = i + 1;
                Complex k0 = Wave.K(waveGuide, m + 1);
                //Complex gamma = Complex.Sqrt(k0 * k0 - (m * Math.PI / waveGuide.A) * (m * Math.PI / waveGuide.A));
                Complex gamma = k0;
                startValues[i + 1] = -gamma * amps[i] * im;
            }
            return startValues;
        }

        public Complex Ey(Complex[] coeff, double x)
        {
            Complex fieldValue = 0;
            for (int i = 0; i < coeff.Length; i++)
            {
                if (i % 2 == 0)
                {
                    fieldValue += coeff[i] * Phi(x, i / 2 + 1);
                }
            }
            return fieldValue;
        }

        public Complex[] Rejection(Complex[,] Ey)
        {
            Complex[] y = new Complex[Ey.GetLength(1)];
            int lastIndex = Ey.GetLength(0) - 1;
            for (int i = 0; i < y.Length; i++)
            {
                y[i] = Ey[lastIndex, i];
            }
            return Rejection(y);
        }

        public Complex[] Rejection(Complex[] y)
        {
            Complex[] rejection = new Complex[y.Length / 2];
            for (int i = 0; i < rejection.Length; i++)
            {
                Complex e = y[2 * i];
                Complex de = y[2 * i + 1];
                Complex k0 = Wave.K(waveGuide, i + 1);
                //Complex gamma = Complex.Sqrt(k0 * k0 - ((i + 1) * Math.PI / waveGuide.A) * ((i + 1) * Math.PI / waveGuide.A));
                Complex gamma = k0;
                rejection[i] = -(de + im * gamma * e) / (de - im * gamma * e);
            }
            return rejection;
        }

        public Complex[] Transmission(Complex[,] Ey, Complex[] rejection)
        {
            Complex[] amps = new Complex[Ey.GetLength(1)];
            int firstIndex = 0;
            for (int i = 0; i < amps.Length; i++)
            {
                amps[i] = Ey[firstIndex, i];
            }

            Complex[] y = new Complex[Ey.GetLength(1)];
            int lastIndex = Ey.GetLength(0) - 1;
            for (int i = 0; i < y.Length; i++)
            {
                y[i] = Ey[lastIndex, i];
            }
            return Transmission(amps, rejection, y);
        }

        public Complex[] Transmission(Complex[] amps, Complex[] rejection, Complex[] y)
        {
            Complex[] transmission = new Complex[rejection.Length];
            for (int i = 0; i < transmission.Length; i++)
            {
                Complex e0 = y[2 * i];
                Complex A0 = e0 / (1 + rejection[i]);
                transmission[i] = amps[i] / A0;
            }
            return transmission;
        }

        public double[] Balanse(Complex[] rejection, Complex[] transmission)
        {
            int length = Math.Min(rejection.Length, transmission.Length);
            double[] balanse = new double[length];
            for (int i = 0; i < balanse.Length; i++)
            {
                balanse[i] = Math.Sqrt(rejection[i].Magnitude * rejection[i].Magnitude + transmission[i].Magnitude * transmission[i].Magnitude);
            }
            return balanse;
        }
    }
}
